#include <iostream>
#include "includes/dg-evolve-calculate-flux.hh"
// #include "includes/basis-function.hh
DGEvolveCalculateFlux::DGEvolveCalculateFlux(){
};
DGEvolveCalculateFlux::~DGEvolveCalculateFlux(){
};
void DGEvolveCalculateFlux::dg_evolve_calculate_flux(SubprojectBase* subproject, double* globalcoord, double* unitouternormal, double pillarlen, int i, int j, int currentdim, int* interscounter, double rdt){
  double unitouternormalsum=0.0;
  double var_i_left[memory::phynumber];
  double var_i_right[memory::phynumber];
  double var_i_middle[memory::phynumber];  
  double var_j_left[memory::phynumber];
  double var_j_right[memory::phynumber];
  double var_j_middle[memory::phynumber];  

  double fequ_i_left[memory::phynumber];
  double fequ_i_right[memory::phynumber];
  double fequ_i_middle[memory::phynumber];  
  double fequ_j_left[memory::phynumber];
  double fequ_j_right[memory::phynumber];
  double fequ_j_middle[memory::phynumber];  

  //double basis_left[memory::dgorder];
  //double basis_right[memory::dgorder];
  //double basis_middle[memory::dgorder];
  double basis_left[8];
  double basis_right[8];
  double basis_middle[8];    
  double flux[memory::phynumber*memory::dgorder];

  double var_fromLeft[memory::phynumber];
  double var_fromRight[memory::phynumber];
  double slope_i[memory::phynumber];
  double slope_j[memory::phynumber];  
  
  // std::cout << "dg_evolve begin" << std::endl;
  for (int tmpi=0; tmpi<gridsetting::dimension; tmpi++)
    unitouternormalsum+=unitouternormal[tmpi];

  basis.getBasisLeft(basis_left);
  basis.getBasisRight(basis_right);
  basis.getBasisMiddle(basis_middle);  
  basis.getVarLeft(i,var_i_left);
  basis.getVarRight(i,var_i_right);
  basis.getVarMiddle(i,var_i_middle);
  
  if (j>=0){
    basis.getVarLeft(j,var_j_left);
    basis.getVarRight(j,var_j_right);
    basis.getVarMiddle(j,var_j_middle);    
  } else {
    // flag = -1.0 or 1.0, for left and right respectively
    subproject->boundary->boundary_function(globalcoord,i,unitouternormalsum,var_j_left,-1.0);
    subproject->boundary->boundary_function(globalcoord,i,unitouternormalsum,var_j_right,1.0);
    subproject->boundary->boundary_function(globalcoord,i,unitouternormalsum,var_j_middle,0.0);    
  }

  // std::cout << "dg_evolve moving-indicator" << std::endl;    

  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i_left,fequ_i_left);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i_right,fequ_i_right);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i_middle,fequ_i_middle);    
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_j_left,fequ_j_left);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_j_right,fequ_j_right);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_j_middle,fequ_j_middle);  

  // std::cout << "begin space discreztion" << std::endl;

  // for second-order algorithem
  for (int tmpi=0; tmpi<memory::phynumber; tmpi++){
    slope_i[tmpi]=memory::limit_slope[tmpi][currentdim][i];
    slope_j[tmpi]=memory::limit_slope[tmpi][currentdim][j];  
    if (unitouternormal[currentdim]>0) {
      var_fromLeft[tmpi]=var_i_middle[tmpi]+slope_i[tmpi]*pillarlen/2.0;
      var_fromRight[tmpi]=var_j_middle[tmpi]-slope_j[tmpi]*pillarlen/2.0;
    } else {
      var_fromLeft[tmpi]=var_j_middle[tmpi]+slope_j[tmpi]*pillarlen/2.0;
      var_fromRight[tmpi]=var_i_middle[tmpi]-slope_i[tmpi]*pillarlen/2.0;
    }
  }

  
  // --- simplest algorithm ------
  switch(setting::spacediscrete){
  case setting::lf:
    for(int tmpj=0; tmpj<memory::dgorder; tmpj++)
      for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
  	if (unitouternormalsum<0){
	  //  	  flux[tmpj*memory::phynumber+tmpi] = (fequ_i_left[tmpi]+fequ_j_right[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_left[tmpi]-var_j_right[tmpi]);
	  //  	  flux[tmpj*memory::phynumber+tmpi] *= (2.0*tmpj+1.0)/2.0*pillarlen*basis_left[tmpj];
  	}
  	if (unitouternormalsum>0){
	  // flux[tmpj*memory::phynumber+tmpi] = (fequ_i_right[tmpi]+fequ_j_left[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_left[tmpi]-var_j_right[tmpi]);
	  //  	  flux[tmpj*memory::phynumber+tmpi] *= (2.0*tmpj+1.0)/2.0*pillarlen*basis_right[tmpj];
  	}
	// flux[tmpj*memory::phynumber+tmpi] = (fequ_i_left[tmpi]+fequ_j_right[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_left[tmpi]-var_j_right[tmpi]);
	  //flux[tmpj*memory::phynumber+tmpi] *= (2.0*tmpj+1.0)/2.0*basis_right[tmpj];


	// flux[tmpj*memory::phynumber+tmpi] = (fequ_i_left[tmpi]+fequ_j_left[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_left[tmpi]-var_j_left[tmpi]);
	flux[tmpj*memory::phynumber+tmpi] = (fequ_i_middle[tmpi]+fequ_j_middle[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_middle[tmpi]-var_j_middle[tmpi]);	
	// flux[tmpj*memory::phynumber+tmpi] = (fequ_i_right[tmpi]+fequ_j_right[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_right[tmpi]-var_j_right[tmpi]);
	// flux[tmpj*memory::phynumber+tmpi] = (fequ_i_left[tmpi]+fequ_j_right[tmpi])/2.0+unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_right[tmpi]-var_j_right[tmpi]);	
      }
    break;
  case setting::upwind:
    break;
  }
  
  // --- high-order algorithm ------  
  switch(setting::spacediscrete){
  case setting::tvdlf:
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      for(int tmpj=0; tmpj<memory::dgorder; tmpj++){	
	flux[tmpj*memory::phynumber+tmpi] = (fequ_i_middle[tmpi]+fequ_j_middle[tmpi])/2.0-0.5*pillarlen/rdt*timerule::cfl/gridsetting::dimension*(var_fromRight[tmpi]-var_fromLeft[tmpi]);
	//	  +unitouternormalsum*pillarlen/(rdt*2.0*gridsetting::dimension)*(var_i_middle[tmpi]-var_j_middle[tmpi]);
      }
    }    
    break;
  case setting::muscl:
    break;
    //  case setting::weno:
    //    break;    
  }

  // caculate update
  if (i==0)
    std::cout << "update with flux" << std::endl;
  for(int tmpj=0; tmpj<memory::dgorder; tmpj++)
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      memory::update[timerule::roundnumber][tmpj*memory::phynumber+tmpi][i] -= 1.0/pillarlen*unitouternormalsum*flux[tmpj*memory::phynumber+tmpi];  

  // std::cout << "update with flux, done" << std::endl;

  if (i==0)
    std::cout << "dg evolve for flux calculation" << std::endl;
};
